#ifndef QUADRIC_H
#define QUADRIC_H


#include "definitions.h"

// This class provides the implementation of the following quadric.
//
//      |  a^2     ab      ac      ad   |
//      |  ab      b^2     bc      bd   |
//      |  ac      bc      c^2     cd   |
//      |  ad      bd      cd      d^2  |
//
// Being the matrix 4x4 symmetric, only the coefficients of the upper triangle are stored
//
class Quadric
{


                public:


                                // Construct quadric from upper triangle of symmetrix 4x4 matrix
                                //
                                Quadric(double _a, double _b, double _c, double _d,
                                                   double _e, double _f, double _g,
                                                              double _h, double _i,
                                                                         double _j) : a(_a), b(_b), c(_c), d(_d),
                                                                                             e(_e), f(_f), g(_g),
                                                                                                    h(_h), i(_i),
                                                                                                           j(_j)
                                {}



                                // Constructor quadric from given plane equation: ax + by + cz + d = 0
                                //
                                Quadric(double _a=0.0, double _b=0.0, double _c=0.0, double _d=0.0) :  a(_a*_a), b(_a*_b),  c(_a*_c),  d(_a*_d),
                                                                                                                 e(_b*_b),  f(_b*_c),  g(_b*_d),
                                                                                                                            h(_c*_c),  i(_c*_d),
                                                                                                                                       j(_d*_d)
                                {}



                                // Construct from point and normal specifying a plane
                                //
                                Quadric(const OpenMesh::Vec3f& normal, const OpenMesh::Vec3f& point)
                                {

                                        *this = Quadric(normal[0], normal[1], normal[2], -dot(normal,point));
                                }



                                // set all matric entries to zero
                                //
                                inline void clear()  { a = b = c = d = e = f = g = h = i = j = 0.0; }



                                // add a given quadric to this quadric
                                //
                                inline Quadric& operator+=( const Quadric& _q )
                                {

                                        a += _q.a;  b += _q.b;  c += _q.c;  d += _q.d;
                                                    e += _q.e;  f += _q.f;  g += _q.g;
                                                                h += _q.h;  i += _q.i;
                                                                            j += _q.j;
                                        return *this;
                                }



                                // multiply the quadric by a scalar
                                //
                                Quadric& operator*=( double _s)
                                {

                                        a *= _s;  b *= _s;  c *= _s;  d *= _s;
                                                  e *= _s;  f *= _s;  g *= _s;
                                                            h *= _s;  i *= _s;
                                                                      j *= _s;
                                        return *this;
                                }



                                // evaluate the quadric Q at position p by computing (p^T * Q * p)
                                //
                                double operator()(const OpenMesh::Vec3f& p) const
                                {

                                        const double x(p[0]), y(p[1]), z(p[2]);

                                        return a*x*x + 2.0*b*x*y + 2.0*c*x*z + 2.0*d*x +
                                                           e*y*y + 2.0*f*y*z + 2.0*g*y +
                                                                       h*z*z + 2.0*i*z +
                                                                                     j;
                                }




                private:

                                double a, b, c, d,
                                          e, f, g,
                                             h, i,
                                                j;
};


#endif // QUADRIC_H
